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Abstract 

Image smoothing and segmentation algorithms are frequently formu- 
lated as optimization problems. Linear and nonlinear (reciprocal) re- 
sistive networks have solutions characterized by an extremum princi- 
ple. Thus, appropriately designed networks can automatically solve 
certain smoothing and segmentation problems in robot vision. This 
paper considers switched linear resistive networks and nonlinear resis- 
tive networks for such tasks. The latter network type is derived from 
the former via an intermediate stochastic formulation, and a new result 
relating the solution sets of the two is given for the "zero termperature" 
limit. We then present simulation studies of several continuation meth- 
ods that can be gracefully implemented in analog VLSI and that seem 
to give "good" results for these non-convex optimization problems. 
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1 Introduction 

One of the most important, yet most difficult, early vision tasks is that of image 
smoothing and segmentation. Smoothing is necessary to remove noise from an input 
image so that reliable processing in subsequent stages is facilitated. However, indis- 
criminate smoothing will blur the entire image, including edges (e.g., corresponding 
to object boundaries) which are necessary for later stages of processing. Many re- 
searchers are currently seeking to develop algorithms that smooth in a piecewise 
manner, respecting edges. There are two main approaches taken — stochastic, [1] - 
[6], and deterministic [7] - [9]. The former relies on such methods as simulated an- 
nealing to accomplish the minimization. The deterministic approach, on the other 
hand, often relies on the application of continuation methods [2], [10] to certain 
nonlinear systems, or in the case of [11], on using a neural network similar to that 
of Tank and II op field [12]. 

Although efficient computation techniques exist for numerically computing the 
solutions to vision problems [13], even the fastest algorithms running on a parallel 
supercomputer (such a,s the Connection Machine ® system 1 [14]) do not approach 
real-time performance. The motivation of this work is to produce solutions to the 
smoothing and segmentation problem that are amenable to analog VLSI network 
implementation, an area that has been explored in [15] - [18]. See also [11], [19], 
[2UJ. 

Section 2 presents the smoothing and segmentation task as a minimization prob- 
lem. Section 3 presents methods for solving the minimization problem and discusses 
network implementations of these methods. Simulation results are provided in Sec- 
tion 4. Finally, conclusions and suggestions for further research are given in Section 
5. 
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2 Image Restoration as a Minimization Problem 

The difficulty with using a linear network for image smoothing is that noise and 
signal are equally smoothed so that edges become blurred. We therefore seek a 
method for segmenting the signal into regions which can be smoothed separately. 
One technique for doing this is to introduce a line process (i.e., a set of binary 
variables) which selectively breaks the smoothness constraint at given locations. 
This method appears widely in the literature, e.g., [1] - [4], [6], [11]. 

For simplicity of notation, all equations in this paper are formulated for the 
one-dimensional case. The results generalize trivially to two dimensions, and the 
simulation results are for the two-dimensional case. 

The smoothing and segmentation problem with the line process can be treated 
as a minimization problem. Let u G $ N be the input image, y e N N be the output 
image, and 1 £ Jc^ -1 be the line process, where the binary line process variable /,• 
assumes the values {0,1} depending on whether the smoothness penalty between 
nodes i and i + 1 is enforced or not. Consider the following cost function: 

«Ai(y>I)= g Wi(y) + 5(y, 1) + //(1)] (1) 

where F, S, and // arc the "fidelity," "smoothness," and "line" penalty terms, 
respectively, i.e., 

A 



^u(y) = *f52(yi-ui) 2 ( 2 ) 

A 7V ~ 1 

S(y,l) = *s52(yi-yi + i)\l-li) (3) 

i'=l 

N-l 

//(I) = XkZh. ( 4 ) 



This formulation assumes that the optimal reconstructed image and edges {y opt ,l 0}t ) 
satisfies as well as possible the generally conflicting requirements of agreeing \vith 
the data u, being smooth between edges, and containing as few edges as possible. 
The parameters A/, A s , and X h determine the weights given to each of these criteria, 
^ The expression (1) can be minimized with respect to y for fixed 1 by differenti- 
ating with respect to each y t and setting the derivatives to zero. This produces the 
following system of equations: 

A/fe - ui) + X s ( yi - t/^XI - /,._,) + \ s ( Vi _ y. +1 )^ __ /.) = 0? (5) 
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*) u(i+1) 



Figure 1: A simple smoothing network with switches. The vertical and horizontal 
resistors have conductances A; and X a , respectively. 

with appropriate modifications at the boundaries i = 1 and / = N . 

The new results in Section 2.1 all follow from the observation that (5), as well 
as the related equation (9) below, describe the solution to certain resistive electrical 
networks. Notice that (5) can be viewed as the Kirchoff's current law (KCL) relation 
at every node of a resistive ladder in which the horizontal resistive elements have 
switches (corresponding to aline process clement) associated with them. A network 
for computing y given 1 is shown in Figure 1. Note that (F + S) is the electrical 
power dissipated in the resistors. Similar networks have appeared in [6] and [1 1]. 
This type of network will be referred to as a resistor- with-switch (or RWS) network. 
For any setting of the switches, the network solution automatically minimizes the 
cost function with respect to y. The difficulty is in minimizing with respect to 1. 

Much has been said in the literature in regard to finding a. global minimum to 
(1) by stochastic and deterministic methods. These techniques are necessary to 
find the minimizing 1 — minimizing with respect to y given 1 only requires the 
solution of a linear system. The deterministic approaches rely on the fact that the 
minimization problem can be recast into one in which the line process variables 
have been eliminated. The latter will be studied here since they appear to lead to 
practical VLSI implementations. 



2.1 Discontinuous Resistive Fuse Elements 

The line process variables can be removed from (1) by straightforward algebraic 
manipulations. In fact, Blake and Zisserman [7] demonstrated that the original 
cost function ./ u (y,l) containing real and boolean variables is intimately related to 
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the following cost function containing only real variables: 



*u(y) 1 1 



N N-l 



(6) 



where 



G(„) i .„!,, {A„„M -/) + A,,/} = { W. H<>/¥. (7) 

tet ' -• [ A/ n otherwise 

The line process is found « posteriori according to: 




\yi - y l+l \ < 
otherwise? 



:s.) 



Note that K u is a non-convex cost function with respect to y. 

Apart from instances in which solutions occur at points where G is not differen- 
tiable, the minimum of K u is to be found among those points where VA' u (y) = 0, 
i.e., 

h(Vi - ui) + g(Vi - yi-i) + g{yi - y l+ i) = 0, (9) 

where g(v) = 1*-G(v). 

Equation (9) can also be viewed as the KCL relation at each node of a. nonlinear 
resistive network with the topology illustrated in Figure 2. The nonlinear resistor 
characteristic, g(v), is that of a linear resistor that reversibly becomes an open 
circuit when the voltage across it exceeds a certain threshold, as shown in Figure 
3. Then, in electrical terms, G is twice the co-content function for this nonlinear 
resistor [16],[21], i.e., G(v) = 2ftg{u)du. We refer to an element of this type as 
a discontinuous resistive fuse and to a network incorporating resistive fuses as an 
RWF network, i.e., a rcsisior-witJi-j'use network. 

For a given cost function, one can construct corresponding RWS and RWF net- 
works. For every solution of an RWF network, there exists a similar solution to 
the corresponding RWS network, but there are switch configurations of an RWS 
network for which there is no corresponding solution in a corresponding RWF net- 
work. The question then arises whether restricting attention to the RWF network 
might cause one to overlook a solution to the RWS network that is in fact a local 
minimum and therefore of potential use in an optimization procedure. The answer 
is no, by the following proposition: 

Proposition 1 Consider the cost function y„(y,l) as specified in (1) for a one- or 
two-dimensional network and the corresponding RWS and RWF networks specified 
by (5) and (9). If the switches are set so that the solution y* for the RWS network 
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Figurc 2: Resistor- Witli- Fuse (UWF) network topology. 




Figure 3: Characteristic of the discontinuous nonlinear resistor known as a "resistive 
fuse." The dotted vertical lines are not part of the constitutive relati 



ion 



2.1 Discontinuous Resistive Fuse Elements ^ 

is not also a. solution to the RWF network, then y* is not a local minimum to ,/„, 
meaning that changing the setting of a single (appropriately selected) switch in the 
RWS network will produce a new solution (with a new value of y) for which the 
value of J u is strictly lower. 

Remark Proposition 1 differs from the result in [7] (pp. 112-113) in that it 
concerns /oca/minima and applies after the network has settled to a new y following 
the closing of the switch. 

In order to complete the proof, we need the following lemma: 

Lemma 1 (Local Measurement Principle) Consider a one-dimensional or two- 
dimensional network of the type shown in Figure 1, in which an arbitrary number 
of switches (> 1) are open. The power dissipated in the network is given by 

N N-i 

P = X ! X>< - m? + K £>,- -y i + x ) 2 0-/,) 
i-i j =1 

= ( F +V- (10) 

For a given switch, let ./>- be the value of P with the switch open, P+ bo the 
value of P with the switch closed, and define the change in the dissipated power as 
^P - P - P~ . Let v oc be the voltage across the switch when it is open and let i NC 
be the current through the switch when it is closed (after the network has settled)" 
Then the increase in dissipation which results from closing the switch is 

&P = v oc i sc >0. (11 ) 

Remark Lemma 1 is a startling result. The local measurement principle states 
that one can measure the (jlobal change in the network cost function due to a switch 
change (after the network settles to a new solution) merely by taking two mea- 
surements at the switch. Both proofs below use circuit theory techniques, but can 
also be carried out, albeit laboriously, by mathematical arguments divorced from a 
network realization, e.g., the proof of Lemma 1 via a rank one perturbation method 
in [22]. Related work appears in [23] and [24]. 

Proof of Lemma 1 Define ?;~ and i~ to be the network branch voltages and 
currents when the switch is open. Define v k and i k to be the network branch 
voltages and branch currents when the switch is closed. Define Av k = v k ~ v~ and 
Aifc = i k - i^. By Tellegen's theorem [25], [26], 



J2 [vk^ik - ik&v k ] = 0. 



2) 



all 
branch 
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Group the terms in (12) according to branch element, and note that 

£ [tfcAt* - tfcAv*] + £ [vkAik - ik&v k ] (13) 



voltage 
sources 



+ X) i v k Ai k ~ h&v k ] + v sw Ai sw - i sw Av sw = 0, 



fixed 
switches 



where the subscript "sw" refers to the switch that is being closed and "fixed 
switches" to all others. To simplify (13), note that Av k = for the voltage sources, 
v k and Av k vanish for closed switches, i k and Ai k vanish for open switches, and for 
the resistors: 

v k Ai k - i k Av k = R k i k Ai k - i k R k Ai k = 0. (14) 

Equation (13) then becomes 

= X l v k^k] + v sw Ai sw - i sw Av sw 

voltage 
sources 



- ]T [ v k&h] + v aw (i aw - i 3W ) - i sw (v sw - v - w ) 

voltage 
sources 

= J2 l V k&h] + V~ w i sw 



voltage 
sources 



(15) 



The summation term in (15) is just the change in power delivered to the network, 
i.e., -AP, and v~ w i sw = v oc i sc . Therefore, 



AP = v oc i sc . 



(16) 



Proof of Proposition 1 Consider any RWS network with any input u, switch 
configuration 1, and corresponding network solution y*, such that y* is not a solution 
of the corresponding RWF network. Then there must exist some resistor-switch 
composite element (element q, say), such that y* is no longer a network solution 
if a resistive fuse is substituted in its place. Make such a substitution and then 
consider the load-line describing the remainder of the linear RWS network as seen 
from this location. The two possible cases marked in Figure 4 are Case X, in which 
switch q was open in the original RWS network, and Case W, in which switch q was 
closed. Note that the area in the first quadrant under the triangle is ±\ h . In Case 
X, closing the switch in the original RWS network would have caused the solution 
to move to the circled point on line A. By Lemma 1 the change would be 

Ju,closed - J\i,open = -[v oc isc — Afc] < 0, (17) 
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voltage 



Figure 4: Load-Line diagram for RWS network and for RWS network with resistive 
fuse substituted for one resistor-switch composite element. The dashed lines A, B, 
and C are possible load-lines representing the behavior of an RWS network as seen 
by one resistor-switch pair. The six marked points indicate possible solutions, de- 
pending on switch position. If a resistive fuse element (with characteristic i = g(v), 
shown with a solid line) is substituted for the resistor-switch pair, the four circled 
solutions remain, while the solutions marked X and W disappear. 

where the inequality follows from the fact that \v oc i sc (the area under the line 
connecting the origin to (v oc ,i sc )) is less than ±\ h (the area under the triangle). 
For Case W, similar reasoning shows J U)0pen - J^ dosed < if opening the switch 
causes the network solution to move from point W to the circled point on line C. 
Thus points X and W in Figure 4 are not local minima of J u . | 

Remark The converse of the proposition is not true. If the network solution lies 
on load-line B, one intersection point or the other will generally have lower cost for 
the RWS network, yet both are valid solutions to the RWF network. 



2.2 Continuous Resistive Fuse Elements 

We will show that for the purposes of numerical optimization and physical (VLSI) 
implementation it is advantageous to replace the discontinuous resistive fuse ele- 
ment in Figure 3 by a controllable element with a single-parameter family of i-v 
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curves, such as those drawn in dotted lines in Figure 7. Elements of this general 
type have also been used in analog VLSI circuits for image enhancement. To the 
best of our knowledge, the first related circuit appeared in [15] and had a mono- 
tone, saturating characteristic of the general form i = atanh(6v). John Harris at 
Caltech has invented the first nonmonotone circuit element of this type, named it 
a resistive fuse, and built image processing networks using it in analog VLSI [16] 
- [18]. More recently, Steve Decker, Hae-Seung Lee, and John Wyatt at MIT have 
developed more compact nonmonotone continuous resistive fuse circuits using fewer 
transistors. 

The behavior of the nonmonotone fuses in the network in Figure 2 is intuitively 
easy to understand. In a smooth region of the image where the input u is nearly 
constant, only the linear portion of the fuse curve near the origin is excited, i.e., 
the fuse acts essentially as a linear smoothing element. But at any point where a 
discontinuity in the input occurs, i.e., where \u { - u^ is sufficiently large, the fuse 
current becomes quite small and little smoothing results. An extremum formulation 
of this behavior is given in [16]. 

A fundamental question is whether any rigorous relationship can be found that 
connects the continuous nonmonotone fuse curves in Figure 7 with the discontinuous 
fuses in Figure 3 or the switches in Figure 2. The surprising answer is yes, due 
to a remarkable result of Geiger and Girosi, based on a stochastic formulation of 
the problem [1]. Section 2.3 gives the necessary background and Section 2.4 gives a 
variant on their approach, based on a formulation in terms of a marginal probability 
distribution function. 

2.3 Stochastic Formulation of the Image Smoothing and Segmen- 
tation Problem 

This section shows that the deterministic minimization problem in (1) - (4) is in 
fact the optimum image reconstruction procedure in a particular probabilistic for- 
mulation of the problem (see Figure 5). In this formulation, the original image 
brightness and discontinuities are modeled as a pair of random vectors (B,D). We 
cannot directly measure (B,D) but have to work instead with a noisy observation 
U of the brightness values alone. 

Remark The notation in this section follows the one classically used in probability 
theory, where the random variables or vectors are denoted by uppercase letters 
(B,D,U) while the values taken by them are denoted by lowercase letters (b, d,u). 
In (18), p(X = x) is denoted by p(x) and p(X = x|Z = z) is denoted by p(x|z) to 
simplify notation. 

For this analysis, it is assumed that the probability distribution for (B,D) is 
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Figure 5: The original image one wishes to reconstruct is modeled as a random 
brightness vector B and a random binary edge discontinuity vector D, where D{ = 1 
if there is an edge between pixels i and (t + 1). The observation U is a version of 
B degraded by observation noise N, and the reconstruction algorithm produces 
an estimate (B,D) of (B,D). In our particular example, the function / is the 
optimization procedure specified in (23) - (25). 



given by 



p(b,d) = p(b|d)p(d) 



(18) 



where the functions S and // are given in (3) and (4), respectively, and where c s 
and c h insure that p(b,d) and p(d) have unit area. The vector D represents a 
finite Bernoulli sequence. More specifically, the components of the binary random 
vector D are independent, identically distributed binary random variables with a 
probability mass function 



so that the joint probability mass function of d is given by 



N-l 



pp=d) = n i 

1=1 



-0Xhd, 



+ e~P x h 



(19) 



(20) 



( 1 + e -^f-^ C ^^ //(d) - 
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which decreases with the total number of discontinuities in the scene. The brightness 
vector B is a Gaussian random vector dependent on D: specifically, the elements 
of the vector taken as a sequence represent, in the one-dimensional case, discrete 
Brownian motion with a uniformly distributed initial value at i = 1 and at the right 
of each discontinuity location. In the absence of a discontinuity between pixels i 
and (i + 1), the standard deviation of B i+1 - B { is a = 1/^/2^7. 
We also assume the observations U are distributed as 

p(u|b) = C/ (^A / )c-^( b ), (21) 

where the function F u is given in (2) and where c f normalizes p(u|b) to unit area. 
Equation (21) is equivalent to assuming that the observations are corrupted by 
additive, independent Gaussian noise, i.e., U { = B x + TV;, where N { is a Gaussian 
random variable with zero mean and standard deviation a = If^jJX], and N is 
independent of (B,D). The variable (3 is actually redundant because scaling is 
equivalent to scaling A,, A s , and \ h . Since increasing (3 will reduce all the variances, 
£ is analogous to temperature in statistical mechanics. 

Remark More precisely, B is Brownian motion and N is Gaussian if the allowed 
values for each B { and each N { are continuous and unconstrained. But in that case 
the description of B is flawed because a uniformly distributed initial value over the 
whole real line is not a normalizable probability distribution. This problem vanishes 
if the allowed brightness values are discrete and finite in number or continuous and 
bounded, e.g., < B x < B max . 

How should we best attempt to reconstruct (B,D) given U? Given both the 
a priori (i.e., prior to a particular observation) probability distribution p(b,d) 
of the image with discontinuities, and the noise-produced conditional distribution 
p(u|b, d) of the observation, Bayes rule [27] gives the a posteriori distribution (i.e., 
the conditional distribution after noisy observation) of the image and discontinuities 
by the formula 

ru j i ^ K u |b,d)p(b,d) 
rtM ' U) = p(u) • ( 22 ) 

^ Maximum a posteriori estimation is a reconstruction technique that chooses 
(b,d) as the value of (b,d) that maximizes p(b,d|u), i.e., 

(b,d) = /(u) = argmaxp(b,d|u) 

b,d ' ' 



— arg max 



>(u|b,d)p(b,d) 



b,d L p(u) 

= argmax[p(u|b)p(b,d)], (23) 
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where the last equality holds because the denominator is independent of (b,d) 
and the observation u depends directly on the brightness levels alone — d is not 
measured but only statistically inferred. Using (18) and (21) in the last line of (23), 
we have 

(b,d) = argmax{c / c.c fc c- /? ^( b )+ 5 ( b ' d )+^( d )]} . ( 24 ) 

Remark There is a useful distinction between the (y,l) notation and the (b,d) 
notation. In the deterministic picture, (y,l) represents algorithm or circuit variables 
over which we are attempting to optimize. A particular network solution (y, 1) does 
describe the circuit behavior but may or may not bear any simple relation to (b, d). 
The stochastic picture adds a new quantity not present in the earlier determinis- 
tic story: the original uncorrupted random image-discontinuity pair (B,D). The 
variables (b,d) always refer to possible values of (B,D) and may or may not relate 
directly to circuit behavior. Without this dual notation the variables would mis- 
leadingly be used to describe both original image-discontinuity pairs and also node 
voltages and switch positions inside an electrical network. 

Substituting our previous deterministic notation (y,l) for (b,d), we recover 

(b,d) = (y op tXpt) = argmin{F u (y) + £(y,l) + //(l)} (25) 

as in (1) - (4). 

In conclusion, the optimization problem in (1) - (4) yields as its solution the 
maximum a posteriori estimator of the original image brightness and discontinuity 
vectors (B,D), assuming the a priori distribution (18) and the observation noise 
model (21). 

2.4 Derivation of the Continuous Fuse from the Resistor-with- 
Switch Network in a Probabilistic Formulation 

In the Bayesian formulation, one first calculated the a posteriori distribution 

p(b, d|u) = C/CsC/l e-^ F u( b )+ 5 ( b ' d )+"( d )], (26) 

and then attempted to maximize it over (b,d). If one wishes to reconstruct only 
the intensities but not the discontinuity locations, it is appropriate to maximize the 
simpler marginal a posteriori distribution p(b\u) over b, where 



A 



p(b|u) = £>(b,d|u), (27) 

dec 



14 
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and C is the hypercube of all 2 N ~ 1 possible binary d-vectors. The form of this 
density will specify the nonlinear continuous fuse characteristics. The sum above 
was first calculated by Geiger and Girosi in [1], and a somewhat more detailed 
derivation is given below. 

Lemma 2 The marginal a posteriori distribution p(b|u) over b is given by: 

p(b|u) = Cie -^ Fu ( b )+ J 2( b )] (28) 

where ci is a normalizing constant, F u (b) is given in (2), and 

N ' x ' 1 + e^ 



1 N ~ x ( 



+ e /3[A h -A,(6,-6, +1 ) 2 ] 



(29) 



The proof of the lemma requires the following fact, which can easily be verified. 
Fact Let a = (a u ... ,a n ) be a vector of n binary variables, a { G {0, 1}, and let A 
be the set of all such vectors. Then for any r£ JJ", 



£«"■' = II(i + «"), 



(30) 



where a • r is the standard inner product. 
Proof of Lemma 2 

Substituting (26) into (27) yields 



Kb|u) = £>(M|u) 



dec 



= c f c s c h J2 e -/3[Fu(b)+S(b,d)+//(d)] 
dec 
The terms being summed in (31) can be decomposed as follows: 



(31) 



exp -(3 



N N-l 

A /E( 6 «-- tt «') 2 + A. E^-W 



N-l 



ex P -^ E d i [ X h ~ Uk - b t+1 ) 2 } .(32) 



i=i 



Using (30), the second exponential term in (32) sums over d G C to 
exp 1-/3 ~jjJ2 Ml + e-fl A *-M 6 -- 6 '+i) 2 ]) 

and further algebraic manipulation shows that 

p(b|u) = c/c.Cfce-^bH^b)^ 



(33) 



(34) 
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where 



N-l 

(35) 



1 / 1 \ 

J 2 (b) = - g In { , + eP[Xh _ Xs{bt _ bt+iy] ) + (N - l)X h . 



Absorbing an additive term into the normalizing constant c f c s c h , the lemma was 
stated in (29) in terms of 

Mb) = J 2 (b) + *Lll [in (l + e^) - 0\ h ] , (36) 

which is constructed so that J 2 (0) = 0. This is a necessary step if we are to later 
interpret J 2 (b) as the co-content function of a set of nonlinear resistors. I 

The marginal distribution in (28) suggests a new cost function 

^(y) = i[^u(y) + J 2 (y)], (37) 

in which the line process variables have been eliminated. The local minima of A'£(y) 
are obtained from the set of points satisfying 

VA'£(y) = 0. (38) 

Taking the i-th component of (38) gives: 

A/(ft - «i) + 9d(y t - y,-_i) + gpiyi - y i+1 ) = 0, (39) 

where 

aM = i + jtv) - m 

Equation (39) can be considered the KCL relation at each node of a nonlinear 
network having vertical linear resistive elements with conductance X f and horizontal 
nonlinear elements with constitutive relation i = 9(3 (v). In this case, A^(y) is the 
total co-content of the network. Notice that as -> oo, we recover the RWF network, 
i.e., A' u (y) -+ A' u (y). Moreover, we have defined a family of /^-dependent resistive 
elements, illustrated in Figure 7, that can be used in continuation methods. 
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3 Solution Methods 

The resistive fuse and marginal distribution approaches produced switch-free nonlin- 
ear networks with identical topologies (see Figure 2) but with different constitutive 
relations for the nonlinear elements. For either network, multiple solutions generally 
exist. On the theoretical side this is a difficulty because we are trying to find the 
global minimum of a specific cost function. On the practical side this is a difficulty 
because the solution that is obtained by a physical network realization will depend 
strongly on such things as parasitic capacitances and other characteristics of the 
network over which we have little control. We therefore seek some modification of 
the network that will allow us to exercise some control over the solution it finds. 
In this section, we apply continuation methods to the nonlinear smoothing and 
segmentation networks. 

3.1 Example — A Special Case 

The simplest special case that nonetheless provides insight into the phenomenon 
of multiple solutions is the response of a one-dimensional network to a step edge 
input, i.e., 



Ui = 



Uhi i < k 



v Uio i>k ' W 

for some k < N. This corresponds to a step of u hi - u lo > between nodes k and 
k + 1 and serves as a model for the simplest two-dimensional edge, i.e., a step that 
extends across the entire network and is parallel to one of the network "axes." 

For the step input described above, the one-dimensional network has a simple 
circuit equivalent, shown in Figure 6. The simplification proceeds as follows. First, 
we assume that the signal is "well-smoothed" on either side of the step so that each 
nonlinear element can be replaced by an equivalent linear resistance whose value is 
the incremental resistance of the nonlinear element about zero volts. The network 
elements on either side of the step are then replaced by their Thevenin equivalents, 
which are combined into a single linear element and voltage source. The simplified 
network will be referred to as the zero-dimensional case. Analysis of the behavior 
of the network to a step input is reduced to solving the KCL equation at one node: 
some insight into the circuit behavior can be gained by using load-line techniques 
(see Figures 6-8). 

This "linear load-line assumption" holds exactly only for the RWS network with 
fixed switch positions and for the marginal distribution network with /3 = 0. For 
the RWF network and for the marginal distribution network with (3 -+ oo, it is 
exact over the limited voltage range in which no new discontinuities are introduced 



18 



3 SOLUTION METHODS 




Figure 6: Thevenin equivalent circuit for one-dimensional nonlinear smoothing and 
segmentation network with step input. 



into y. Otherwise, it is only an approximation and its applicability to other cases 
of interest must be individually determined. 

3.2 Continuation Methods 

We seek a modification to the networks so that the solution will be repeatable and 
also be visually and quantitatively "good." One technique that works well within 
the context of smoothing and segmentation is to apply a continuation method to 
the network [2], [10]. 

A continuation (sometimes called "deterministic annealing") can be realized in 
network form by the simultaneous application of a given homotopy (continuous 
deformation) to some or all of the circuit elements. Two types of continuations 
are particularly appropriate for our class of nonlinear networks. Assume we have a 
network with horizontal nonlinear resistors whose constitutive relation is described 
by i = g(v), and vertical linear resistors with conductance X f . Consider the following 
two homotopies for the horizontal and vertical elements, respectively: 

CH: Replace g with 0<"), p e [a, 6], such that g( a ) constrains the network to have a 
unique solution and that g^ = #; 

CV: Replace A, with A^ p) , p e [a, 6], such that A^ a) constrains the network to have 
a unique solution and that A^ = A/. 

Note that CH and CV define where the homotopies are applied in the network to 
produce a continuation; we are still free to decide the specific form of the homotopy. 



3.3 P-Continuation 



3.3 /^-Continuation 
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Blake and Zisserman suggest a CH continuation method — the so-called "graduated 
non-convexity" algorithm, or GNC [7]. There are some apparent weaknesses to 
using the GNC algorithm in network form, however. First, there is no reason 
to expect that, for an arbitrary image, the specific continuation used by GNC 
will produce the global minimum or that it will even produce a "good" minimum. 
Second, the nonlinear resistive element in a network realization of GNC will have 
a discontinuous first order derivative which can cause convergence difficulties in 
numerical simulation. 

On the other hand, the marginal distribution derivation of our nonlinear network 
provides a natural homotopy for realization of the CH continuation. For j3 = 0, the 
network with elements described by (40) is linear, whereas for j3 -► oo, the elements 
become identical to those in Figure 3 and will (locally) solve our minimization 
problem. This suggests using p directly as the continuation parameter for a CH 
continuation for solving (39) and hence (9). Furthermore, because of the way this 
continuation was derived, one might expect that it would do a good job of seeking 
the global cost minimum. 

Some insight into the behavior of this type of network can be gained by examin- 
ing the zero-dimensional case. Figure 7a shows the marginal distribution nonlinear 
resistor characteristic for various values of /?, along with two load-lines representing 
two different values of the input. As p is taken from to oo, the solution will follow 
the continuous path represented by the intersection of the resistor curve and the 
load-line. In this example, the smaller step will be smoothed, and the larger step 
will be segmented. 

Interestingly, discontinuous behavior can occur with this type of continuation, 
as is shown in Figure 7b. In this example, the initial solution point will be the 
intersection of the load-line and the marginal distribution resistor characteristic for 
P = 0. As P is increased, the "hump" of the nonlinear resistor curve will at one 
point pass completely beneath the load-line, at which point the solution will jump 
from being a smoothing solution to being a segmenting solution. 

3.4 Ay-Continuation 

The CV continuation can be realized in a straightforward manner by varying the 
vertical resistors in the network. In particular, we begin with the resistors having 
infinite (or sufficiently large) conductance so that the network has only one solution, 
namely y = u (or, for large conductance, y « u). Then, we continuously decrease 
the value of the conductance to Xj. 
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Figure 7: Approximate load-line plot for marginal distribution network with 
/^-continuation. In part (a), the solid lines with negative slope represent the 
load-lines for two different input values (0.5 V and 2.5V). The nonlinear resis- 
tor is shown for various values of (3. For = 0, the nonlinear resistor acts as a 
linear resistor. As /? -» oo, the nonlinear resistor characteristic becomes that of the 
discontinuous resistive fuse. In part (b), the solution exhibits a discontinuous jump, 
from a smoothing solution to a segmenting solution. 



3.4 Xf -Continuation 91 

Examination of the zero-dimensional case provides some insight into the behav- 
ior of this type of network. Figure 8a shows the marginal distribution nonlinear 
resistor characteristic for large /?, along with two sequences of load-lines represent- 
ing two different values of the input. In this example, the solution for the larger 
input will remain at the initial intersection point of the load-line and the resistor 
curve as \f is taken from A^ a) = A to A^ 6) = X f . On the other hand, the solution 
for the smaller input will follow the continuous path represented by the intersection 
of the resistor curve and the load-line. Hence, the larger step will be segmented and 
the smaller step will be smoothed. 

Discontinuous behavior can also occur with this type of continuation, when the 
continuation is used with non-linear resistors of finite /?, as is shown in Figure 8b. 
In this example, the initial solution point will be a segmenting solution in the lower 
right-hand corner of the figure. As A/ is decreased, the load-line will at some point 
pass completely beneath the nonlinear resistor characteristic, at which point the 
solution will jump from being a segmenting solution to being a smoothing solution. 
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Figure 8: Approximate load-line plot for Ay-continuation. In part (a), the nonlinear 
resistor characteristic g p (v) is shown for large along with two sets of load-lines, 
each set for a different value of the input (the load-lines intersect the g p (v) = 
line at the value of the input voltage: 0.5 V and 2.5 V). As X f is decreased, the 
load-lines rotate counter-clockwise. In part (b), the non-linear resistor characteristic 
is shown for finite p. In this case, the solution exhibits a discontinuous jump, from 
a segmenting solution to a smoothing solution, as X f is decreased. 
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Figure 9: TwoDimensional network topology. 



4 Numerical Experiments 

In order to quantitatively and qualitatively demonstrate the behavior of the /3- 
continuation and Ay-continuation networks, the results of several numerical experi- 
ments are presented. The experiments were all conducted with two-dimensional net- 
works, the topology of which is shown in Figure 9. The experiments were conducted 
using serial and parallel versions of a special purpose circuit simulator developed 
specifically for vision circuits [28] - [30]. 

The continuations were simulated by performing dynamic simulations of the 
networks. In order to add dynamics to the networks, a small parasitic capacitance 
to ground was added at each node such that the time constants of the network were 
much faster than the rate at which the circuit elements were varied to perform the 
continuation. Dynamic simulation of the networks in this way has several advan- 
tages. First, the presence of parasitic capacitances is somewhat more physical and 
will allow the system to perform a gradient descent which will thereby guarantee 
that the network does not settle on a solution which statically satisfies KCL but is 
actually a local maximum of the network cost function [16]. Second, the dynamics 
will insure that the network behavior is well-defined at points where solutions in the 
static case would disappear, as in Figures 7b and 8b. (Our experience has been that 
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(a) 



(b) 



Figure 10: (a) Original image, (b) Original image corrupted with noise. The noisy 
image was used as input for ail experiments in Section 4.1. 

discontinuous circuit behavior is much more common in the Ay-continuation network 
than in the /^-continuation network, which causes simulation of a Ay-continuation 
network to take much more time.) 



4.1 Experiments with a Synthetic Image 

A series of seven experiments conducted on a 16 X 16 circuit grid with a synthetic 
input image was conducted. Figure 10a shows the 16 x 16 synthetic image used for 
the experiments. The small step is 1 V in height and the large step is 3 V. The 
original image was then corrupted by the addition of 0.5 V of uniformly distributed 
noise and is shown in Figure 10b. The noisy signal was used as input for this series 
of experiments. 

For each experiment, a cost function was determined and the corresponding /?- 
continuation and Ay-continuation networks constructed. Then, the networks were 
each simulated using the input image shown in Figure 10b. For each experiment, 
the value of A s was fixed at 1.0 x 10~ 3 and the value of X h was changed. For the' 
/^-continuation, the value of A/ was fixed at lxlO" 4 and the value of (3 was increased 
from to 20/ A^. For the Ay-continuation, the value of/? was set to 20/A^ and the 
value of Ay was varied from 1 to lxl0~ 4 . Thus, for each experiment, the final states 
of the /^-continuation and Ay-continuation networks were the same. 



4.1 Experiments with a Synthetic Image 
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Expt 


A* 


Cost 
/?-cont. Af-cont. 


1 


l.OxlO" 3 


1.775 XlO"' 2 


1.770X10" 2 


2 


5.0 xlO~ 4 


9.699 xlO~ 3 


1.254X10- 2 


3 


1.0 xlO~ 4 


3.299 xlO~ 3 


2.940 XlO" 3 


4 


5.0 xlO" 5 


1.740 xlO~ 3 


1.740 xlO" 3 


5 


l.OxlO" 5 


7.800 XlO" 4 


2.641 xlO~ 3 


6 


5.0 xlO~ 6 


6.600 xlO" 4 


1.650 XlO" 3 


7 


l.OxlO" 6 


5.518 xlO" 4 


4.246 XlO" 4 



Table 1: Experimental results showing the values of the cost function of the solutions 
produced by the (3- and Ay-continuation networks for different values of \ h . 



Solutions obtained by the two nonlinear networks were compared as follows: 

1. Given a cost function, construct the corresponding nonlinear networks, and 
in addition, construct a corresponding RWS network; 

2. Provide each network with the same input and allow each network to attain 
its solution; 

3. For each nonlinear network, transfer the line process solution obtained to the 
RWS network by setting the switches according to equation (8); 

4. Allow the RWS to attain its voltage solution and compute the resulting cost 
— it is this cost that is used for comparison. 

The results of the seven experiments are shown in Table 1. For the particular 
values of parameters used, each network computed a lower cost in roughly half the 
experiments. This set of experiments was actually taken from a larger set of 49. Of 
those, the 0- continuation found the lower cost 35 times, the Ay-continuation found 
the lower cost eight times, and there were six ties. Thus, in these experiments, the 
/^-continuation performs its task of minimizing the cost function (1) extremely well. 

If the cost function were the last word on image smoothing and segmentation, 
we could immediately recommend the /^-continuation. However, remember that the 
ultimate goal for a smoothing and segmentation network is essentially to recover 
an original image minus any noise, and the cost function was introduced to give 
us a quantitative means for doing this. Now consider Figure 11, which corresponds 
qualitatively to the solutions produced by the two nonlinear networks in experiments 
2 and 3. Note that whereas 11a is the qualitatively correct solution, it corresponds to 
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(a) 



(b) 



Figure 11: (a) Network solution produced by Ay-continuation in experiments 2 and 
3. (b) Network solution produced by /^-continuation in experiments 2 and 3. 



the higher cost in experiment two, but corresponds to the lower cost in experiment 
three. 

Naturally, this calls into question the entire cost function methodology used 
for smoothing and segmentation. The difficulty arises because our efforts were 
concentrated only on finding an optimal solution rather than the larger issue of 
determining the best cost function and parameter values. See however [31]. 

4.2 Experiments with a Real Image 

The networks were then tested with a real image. Figure 12 shows the 256 X 256 
input image — a portion of the San Francisco sky line. The output images shown in 
Figures 13 - 16 were produced using a recently developed circuit simulation program 
on the Connection Machine. 

Figure 13 shows the output produced by the /3- continuation with fixed parameter 
values X s = lxlO" 3 , A/, = 2xl(T 5 , and A, = lxlO" 4 . Figure 13a shows the output of 
the network at the beginning of the continuation when (3 = 0; Figure 13b shows the 
output of the network at an intermediate point of the continuation when /3 = 5xl0 3 ; 
Figure 13c shows the output of the network at the end of the continuation when 
= 5xl0 5 . Figure 14 shows the output produced by the /^-continuation with 
fixed parameter values A s = 1 x 10" 3 , X h = lxlO" 5 , and Ay = 3xl0" 5 . Figure 
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14a shows the output of the network at the beginning of the continuation when 
(3 = 0; Figure 14b shows the output of the network at an intermediate point of the 
continuation when (3 = 2x 10 4 ; Figure 14c shows the output of the network at the 
end of the continuation when /3 = 1 x 10 6 . Figure 15 shows the output produced 
by the A/-continuation with fixed parameter values X s = lxlO -3 , A^ = 1 x 10~ 5 , 
and p = lx 10 6 . Figure 15a shows the output of the network at the beginning of 
the continuation when A/ = 1; Figure 15b shows the output of the network at an 
intermediate point of the continuation when Ay = 1 x 10~ 3 ; Figure 15c shows the 
output of the network at the end of the continuation when Xj = 3xl0~ 5 . Note that 
the final parameter values of this network are identical to those for the network 
of Figure 14. Figure 16 shows the output produced by the A/-continuation with 
parameter values \ s = lxl0~ 3 , A^ = 2xl0~ 5 , j3 = 5xl0 4 . Figure 16a shows the 
output of the network at the beginning of the continuation when A/ = 1; Figure 
16b shows the output of the network at an intermediate point of the continuation 
when Xf = 5x 10~ 4 ; Figure 16c shows the output of the network at the end of the 
continuation when Xj = lxlO -6 . 

Discussion 

As can be seen from the experiments with the real image, not only does the se- 
lection of parameter values affect the behavior of the networks, but the continuation 
used also has a profound effect on the network behavior. The differences between 
cost functions for a particular continuation can be seen by comparing Figures 13c 
and 14c, and by comparing Figures 15c and 16c. The differences between continu- 
ations methods for a given cost function can be seen by comparing Figures 14c and 
15c. 

One can understand the differences in the continuation methods quite readily. 
At the beginning of the /^-continuation, the output of the /^-continuation network is 
rather smooth, since initially the network is equivalent to a linear resistive network 
(see Figures 13a and 14a). The edges are then added during the course of the 
continuation (see Figures 13b and 14b). This is a difficulty because without any 
initial edge information, some of the edges might be misplaced or even completely 
lost. Notice that in Figures 13c and 14c, edges tend to line up along the network 
axes. 

On the other hand, the initial output image of the A/-continuation network 
is very close (or identical) to the input one (see Figures 15a and 16a). All the 
edges are initially present and only the spurious edges are smoothed during the 
course of the continuation. Since all the edge information is present at the start of 
the continuation, one would expect that the Ay-continuation would more properly 
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Figure 12: 256 x 256 image of the San Francisco sky line. 



4.2 Experiments with a Real Image 
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Figure 13: (a) Output produced by /^-continuation network. Here, the parameter 
values are X s = lxlO" 3 , A^ = 2xl0" 5 , A, = lxlO" 4 , and = 0, 5xl0 3 , and 5xl0 5 
for Figures 13a, 13b, and 13c, respectively. 
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Figure 13: (b) Output produced by /^-continuation network. Here, the parameter 
values are A s = lxlO" 3 , X h = 2xl0~ 5 , A, = lxlO" 4 , and (3 = 0, 5xl0 3 , and 5xl0 5 
for Figures 13a, 13b, and 13c, respectively. 
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Figure 13: (c) Output produced by /^-continuation network. Here, the parameter 
values are X s = lxlG" 3 , X h = 2xl0" 5 , Ay = lxlO" 4 , and = 0, 5xl0 3 , and 5xl0 5 
for Figures 13a, 13b, and 13c, respectively. 
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Figure 14: (a) Output produced by /^-continuation network with smaller fidelity 
and line penalty weights and larger final (3 value than for the network in Figure 13. 
Here, the parameter values are A 5 = 1 x 10" 3 , A,, = lxlO" 5 , X f = 3xl0" 5 , and 
P = 0, 2 x 10 4 , and 1 x 10 6 for Figures 14a, 14b, and 14c, respectively. 
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Figure 14: (b) Output produced by ^-continuation network with smaller fidelity 
and line penalty weights and larger final /? value than for the network in Figure 13. 
Here, the parameter values are A 5 = lxlO" 3 , A^ = lxlO" 5 , X f = 3xl0~ 5 , and 
(3 = 0, 2 x 10 4 , and 1 x 10 6 for Figures 14a, 14b, and 14c, respectively. 
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Figure 14: (c) Output produced by /^-continuation network with smaller fidelity 
and line penalty weights and larger final (3 value than for the network in Figure 13. 
Here, the parameter values are X a = 1 x 10" 3 , X h = 1 x 10~ 5 , Xj = 3 x 10" 5 , and 
fi = 0, 2 x 10\ and 1 x 10 6 for Figures 14a, 14b, and 14c, respectively. 
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Figure 15: (a) Output produced by A/- continuation network. Here, the parameter 
values are A, = lxlO" 3 , A A = lxlO" 5 , = lxlO 6 , and X f = 1, A, = lxlO" 3 , 
and X f = 3x10 5 for Figures 15a, 15b, and 15c, respectively. Note that the final 
parameter values of this network are identical to those for the network of Figure 14, 
but that the output image is much closer to the input image shown in Figure 12. ' 
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Figure 15: (b) Output produced by Ay-continuation network. Here, the parameter 
values are A 5 =^lxl()- 3 , X h = lxlO" 5 , /3 = lxlO 6 , and X f = 1, \ f = lxHT 3 , 
and A/ = 3X10 5 for Figures 15a, 15b, and 15c, respectively. Note that the final 
parameter values of this network are identical to those for the network of Figure 14, 
but that the output image is much closer to the input image shown in Figure 12. ' 
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Figure 15: (c) Output produced by X f - continuation network. Here, the parameter 
values are A 5 = lxlO" 3 , X h = lxlO" 5 , = lxlO 6 , and A, = 1, X f = 1 x 10" 3 , 
and X f = 3x10 5 for Figures 15a, 15b, and 15c, respectively. Note that the final 
parameter values of this network are identical to those for the network of Figure 14, 
but that the output image is much closer to the input image shown in Figure 12. ' 
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Figure 16: (a) Output produced by A/-continuation network. Here, the parameter 
values are X s = lxlO" 3 , X h = 2xl0" 5 , = 5xl0 4 , and A, = 1, Xj = 5xl0" 4 , and 
X f = lxlO" 6 for Figures 16a, 16b, and 16c, respectively. 
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Figure 16: (b) Output produced by Ay-continuation network. Here, the parameter 
values are A 5 = lxlO" 3 , X h = 2xl0" 5 , fi = 5xl0 4 , and A, = 1, A, = 5xl(T 4 , and 
A/ = 1x10 6 for Figures 16a, 16b, and 16c, respectively. 
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Figure 16: (c) Output produced by A/-continuation network. Here, the parameter 
values are A,= lxlO" 3 , A, = 2xl0~ 5 , = 5xl0 4 , and A, = 1, A/ = 5xl0" 4 , and 
A/ = 1x10 6 for Figures 16a, 16b, and 16c, respectively. 
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locate and preserve edges. This observation is borne out in Figures 15c and 16c; 
the edges are much more well preserved than those in Figures 13c and 14c. 

Finally, since the Ay-continuation only requires that a linear resistance be var- 
ied, its VLSI implementation should be much more compact than that of the j3- 
continuation, which would require varying the characteristics of a nonlinear resistor. 

4.3 Behavior of the Ay-Continuation 

There were some interesting properties exhibited by networks constructed with the 
elements described in (40) having a fixed f3 < oo. For such a network, it can be 
shown that there exist a X min > and a \ max < oo such that for Ay > X max and for 
Ay < X min , the network has a unique solution. In fact, for X f > A max , the output 
will essentially match the input (i.e., y « u), whereas for X f < A mm , the output will 
contain no edges. Consider the network behavior as a function of X f as X f is varied 
continuously from X max to X min . The initial solution of the network will closely 
match the input. Then, as X f is decreased, edges will begin to disappear, first the 
smaller, then the larger, until all the edges are gone. In other words, X f acts as 
a scale-space parameter. This has important practical applications. The dynamic 
network of Perona and Malik [9] has the property that time acts as a scale-space 
parameter. In contrast, we can exercise direct control over the scale-space parameter 
in the Ay-continuation network. See also [2]. Moreover, this behavior is somewhat 
reminiscent of the more successful methods used for hierarchical multiscale image 
representation [32]. 

Under some mild assumptions, it can be shown that the particular solution path 
with the endpoints described above is continuous, connected, and can be numerically 
traced out in U N+1 (Ay X 3l N space) using an arc-length continuation [33]. In 
such a case, any particular value of X f would correspond to an JV-dimensional 
hyperplane parallel to 3l N given by x N+1 = X f and network solutions for this Ay- 
would be intersections of the solution path and that hyperplane. An interesting 
question now arises: why can't we just trace out the path in & N+1 , determine 
the solutions, and sort them by cost to find the global minimum? The answer, 
unfortunately, is that there can exist solution loops that are disjoint from the main 
solution path, meaning that the path traced out from the starting point of large 
Ay will not necessarily contain all the solutions at any given value of Ay. The 
solution loops can occur in as small an example as a three pixel circuit (see Figure 
17) and we offer as proof the experimental evidence in Figure 18. To produce 
the paths, Ay was parameterized as a function of parameter t according to Ay = 
(1 ~ O^iarge + A small , and the plots were made in t x 3c 2 space. Notice that, as 
predicted, there is one solution path with endpoints {t = 0,Vi = V in = 2.5, V 2 = 0} 
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Figure 17: Three pixel example circuit. For the result shown in Figure 18, V tn = 



2.5 V, X s = 0.19, A^ = l.Ox 10- 4 , (3 = 40, A, = (l - i)A large + A s 
Asmaii = l.OxlO -12 , and / was varied from to 1 



mall i 'Marge 



= 0.1, 



and {t = l,Fi = 0,V 2 = 0}, corresponding to an edge across the first nonlinear 
resistor. In addition, there is a closed solution loop centered (approximately) at 
{t = 0.9, Vi = 0.2, V 2 - 1.2}, corresponding to a "misplaced" edge, i.e., an edge 
across the second nonlinear resistor. In general, it can be shown that for a one- 
dimensional RWS or RWF network with a single step input, solutions corresponding 
to misplaced edges always have higher cost than solutions corresponding to correctly 
placed edges. 
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Figure 18: Solution path in t x ft 2 space for the three pixel example circuit demon- 
strating the existence of a disconnected solution loop. 
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5 Conclusion 

In this paper, we developed and compared a series of nonlinear networks for image 
smoothing and segmentation. The results of several experiments indicate that the 
typical cost (or "energy") function minimization formulation of the smoothing and 
segmentation problem does not necessarily capture the essence of the task. For 
the specific parameter values we used, the A/-continuation network performed ex- 
tremely well even though it did not always find the solution with minimum cost. 
The A/-continuation network has several implementation advantages over the (3- 
continuation network. First, in certain cases, it seems to perform the smoothing 
and segmentation task in a more visually correct fashion. Second, X f can be used 
as a scale-space parameter. Finally, since the A/-continuation only requires that a 
linear resistance be varied, its VLSI implementation should be much more compact 
than that of the /^-continuation. 

Several open questions remain. Primary among these is the need for a com- 
prehensive characterization of the natural behavior of these networks. By "natural 
behavior" we mean a set of quantitative empirical statements that relate the behav- 
ior of the network, given certain canonical edge configurations, to the cost function 
parameters. Furthermore, it is important to know how the networks behave in the 
presence of varying amounts and types of noise. Finally, Tom Richardson has devel- 
oped an alternate formulation of the smoothing and segmentation problem based on 
a rigorous analysis of the continuous case [31]. This leads to a more complex circuit 
interpretation that might offer better performance than the methods investigated 
here. Since efficient simulation tools on the Connection Machine are now available, 
it is hoped that some of these questions can be addressed in the near future. 
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